home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / obsolete / spearman.pro < prev    next >
Encoding:
Text File  |  1997-07-08  |  5.5 KB  |  222 lines
  1. ; $Id: spearman.pro,v 1.2 1997/01/15 04:02:19 ali Exp $
  2. ;
  3. ;  Copyright (c) 1991-1997, Research Systems Inc.  All rights
  4. ;  reserved. Unauthorized reproduction prohibited.
  5.  
  6.  
  7.  pro  compute_ave, Data,SortT,R, count
  8.  ; returns a ranking of Data with ties averaged
  9.  ; in SorT. On input sortT = sort(data). On output,
  10.  ; sortt(i) = the ranking of the value in data(i).
  11.  
  12.      temp = Data(SortT)
  13.      temp1 = sortT
  14.      for i = 0,R-1 DO temp1(sortT(i)) = i +1
  15.      sortT = temp1
  16.  
  17.  
  18.    Repeats =  where (temp EQ shift(temp,-1),count) 
  19.                           ; look for repetitions to average
  20.       
  21.       if count NE 0 THEN BEGIN
  22.         here = 0
  23.         
  24.         while ( here LT count) DO BEGIN
  25.  
  26.             first =temp(Repeats(here))
  27.             set1= where( Data EQ first, count2)
  28.             ave = repeats(here) + (count2 - 1.)/2.0
  29.             SortT(set1) = ave 
  30.             here = here + count2 -1
  31.             
  32.        ENDWHILE
  33.       ENDIF
  34. return
  35. end
  36.  
  37.  
  38. pro spearman, Data1, CorMatrix,Names=Names,List_Name=Ln,   $
  39.             Missing=M ,NoPrint = NP
  40. ;+ 
  41. ; NAME:
  42. ;    SPEARMAN
  43. ;
  44. ; PURPOSE:
  45. ;    To compute the correlation matrix where Spearman's correlation 
  46. ;    coefficient is computed between pairs of ranked variables. No 
  47. ;    assumptions are made about the probabilty distribution of the 
  48. ;    variable values.
  49. ;
  50. ; CATEGORY:
  51. ;    Statistics.
  52. ;
  53. ; CALLING SEQUENCE:
  54. ;    SPEARMAN, Data, CorMatrix
  55. ;
  56. ; INPUTS:
  57. ;    Data:    Two-dimensional array. Data(i,j) = the jth value of the ith
  58. ;        variable.
  59. ;
  60. ; KEYWORDS:
  61. ;    NAMES:    Vector of user supplied names for the variables to be used in
  62. ;        the output.
  63. ;
  64. ;    LIST_NAME:    Name of output file.  Default is to the screen.
  65. ;
  66. ;      MISSING:    Value used as a place holder for missing data.  Pairwise 
  67. ;        handling of missing data.
  68. ;
  69. ;     NOPRINT:    A flag, if set, to suppress output to the screen.
  70. ;           
  71. ; OUTPUT:
  72. ;       Correlation matrix written to the screen
  73. ;
  74. ; OPTIONAL OUTPUT PARAMETERS:
  75. ;    CORMATRIX:    Two-dimensional correlation array.  CorMatrix(i,j) = spearman
  76. ;        correlation coeffiecient computed for variables i and j.
  77. ;                          
  78. ; RESTRICTIONS:
  79. ;    None.
  80. ;
  81. ; COMMON BLOCKS:
  82. ;    None.
  83. ;
  84. ; PROCEDURE:
  85. ;    Variables are pairwised ranked with ties averaged.  Correlation 
  86. ;    coefficient is computed using the rankings.
  87. ;-            
  88.  
  89.  
  90.  
  91.  
  92. On_Error,2
  93. Data = Data1
  94. SD= size(Data)
  95. MissNo = 0
  96.  
  97. if( N_Elements( Ln) NE 0) THEN openw,unit,/get,Ln else unit=-1
  98.  
  99. if ( SD(0) NE 2) THEN BEGIN
  100.    printf,unit, 'spearman- Data array has wrong dimension'
  101.    goto, DONE
  102. ENDIF
  103.  
  104.  
  105. C=SD(1)
  106. R= SD(2)
  107.  
  108. CorMatrix = Fltarr(C,C) +1
  109. SortMatrix  = Data
  110. if C mod 2 NE 0 THEN C1 = C+1 else C1 = C 
  111.  
  112. here = where( Data GE 1.0e30 , num)
  113. if num NE 0 THEN BEGIN
  114.   print,unit, ' spearman - data values too large = inf'
  115.   goto,done
  116. ENDIF
  117.  
  118.  
  119. if( N_Elements(M) NE 0) THEN BEGIN         ; put missing data 
  120.                                            ; where it wont
  121.                                            ; affect sort
  122.    here = where(Data EQ M,missNo)
  123.    if (missNo NE 0) THEN Data(here) = 1.0e30
  124. ENDIF ELSE missNo = 0
  125.  
  126.  
  127.  
  128. TN = 0                              ;TN = total number of ties
  129. for i = 0,C-1 Do BEGIN               ;Rank each population
  130.  
  131.    Sort11 =float( sort(Data(i,*)))
  132.    compute_ave, Data(i,*), Sort11, R, TieNo      ; handle ties
  133.    TN = Tn + TieNO
  134.    SortMatrix(i,*) = sort11
  135.  ENDFOR
  136.  
  137. if TieNo EQ 0 and MissNo EQ 0 THEN BEGIN    ;handle simple case
  138.    CorMatrix = 1 - (4*R +2.)/(R - 1.)   $
  139.               + (12.0/(R * (R^2 - 1.)))    $
  140.                          * (SortMatrix # transpose(SortMatrix))
  141.  
  142. ENDIF ELSE   $
  143.   
  144.   for i = 1,C1/2 DO BEGIN           ;cycle through populations
  145.       for j = 0,C-1 DO BEGIN
  146.        Pop = R
  147.        k = (j + (C-i)) mod C
  148.  
  149.        Sort11 = SortMatrix(j,*)
  150.        Sort22 = SortMatrix(k,*)    
  151.  
  152.        if missNo NE 0 THEN  BEGIN    
  153.           here = where ( Data(j,*) EQ 1.0e30,count1)
  154.           if count1 NE 0 THEN BEGIN
  155.              Pop = Pop - count1
  156.              tot_miss = here
  157.              Sort2 = Data(k,*)
  158.              Sort2(here) = 1.0e30
  159.              Sort22 = sort(Sort2)
  160.              compute_ave,Sort2 ,Sort22,R
  161.            ENDIF
  162.  
  163.           here = where ( Data(k,*) EQ 1.0e30 and Data(j,*) NE $
  164.                  1.e30, count)
  165.           if count NE 0 THEN BEGIN
  166.              Pop = Pop - count
  167.              if count1  EQ 0 then   $
  168.                tot_miss = here else  $
  169.                 tot_miss = [tot_miss,here]
  170.              Sort1 = Data(j,*)
  171.              Sort1(here) = 1.0e30
  172.              Sort11 = sort(Sort1)
  173.              compute_ave, Sort1, Sort11, R
  174.            END
  175.         if (Pop NE R) THEN BEGIN
  176.            Sort11(tot_miss) = 0
  177.            Sort22 (tot_miss) = 0
  178.         ENDIF
  179.      ENDIF      
  180.        temp = Pop * Total(Sort11 ^2) - (Total(Sort11))^2 
  181.        temp1 =Pop * Total(Sort22 ^2) - (Total(Sort22))^2
  182.        temp = sqrt(temp * temp1)
  183.        CorMatrix(j,k) = (Pop*Total(Sort11*Sort22) -     $
  184.                         Total(Sort11) *Total(Sort22)) /temp
  185.        CorMatrix(k,j) = CorMatrix(j,k)
  186.     ENDFOR
  187.   ENDFOR
  188.  
  189.  
  190.  
  191.  
  192.   
  193.  
  194.  
  195.    
  196.  
  197.  SN =Size(Names)
  198.  if (SN(1) EQ 0) THEN BEGIN
  199.      I = INDGEN(C)
  200.      Names=['Var'+StrTrim(I,2)]  
  201.  ENDIF ELSE                  $
  202.    if ( SN(1) LT C) THEN BEGIN
  203.      I = Indgen(C)
  204.      printf,unit,'spearman- missing names'
  205.      Names=[Names, 'var'+StrTrim(I(SN(1):C-1),2)]
  206.   ENDIF
  207.  
  208.  
  209.  
  210.   if(N_Elements(NP) EQ 0) THEN BEGIN
  211.     printf,unit, " Correlation Matrix"
  212.     printf,unit, " "
  213.     printf,unit, format ='(10X,16(A8,2x))', Names
  214.    for i= 0,C-1 do                       $
  215.      printf,unit, format='(A8,2X,16(f8.5,2X))',Names(i), CorMatrix(*,i)
  216.   ENDIF
  217.  
  218. DONE:
  219.    if ( unit NE -1) THEN Free_Lun,unit
  220.     RETURN
  221.     END
  222.